home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / public / SciAn / src / ScianVisIso.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  48KB  |  1,933 lines
  1. /*ScianVisIso.c
  2.   Routines for the isosurface visualization object
  3.   Eric Pepke
  4.   May 25, 1991
  5. */
  6.  
  7. #include "Scian.h"
  8. #include "ScianTypes.h"
  9. #include "ScianArrays.h"
  10. #include "ScianWindows.h"
  11. #include "ScianTextBoxes.h"
  12. #include "ScianObjWindows.h"
  13. #include "ScianIcons.h"
  14. #include "ScianColors.h"
  15. #include "ScianControls.h"
  16. #include "ScianButtons.h"
  17. #include "ScianTitleBoxes.h"
  18. #include "ScianLists.h"
  19. #include "ScianSpaces.h"
  20. #include "ScianSliders.h"
  21. #include "ScianIDs.h"
  22. #include "ScianDatasets.h"
  23. #include "ScianErrors.h"
  24. #include "ScianVisObjects.h"
  25. #include "ScianVisWindows.h"
  26. #include "ScianStyle.h"
  27. #include "ScianPictures.h"
  28. #include "ScianDepend.h"
  29. #include "ScianTimers.h"
  30. #include "ScianScales.h"
  31. #include "ScianTemplates.h"
  32. #include "ScianTemplateHelper.h"
  33. #include "ScianEvents.h"
  34. #include "ScianScripts.h"
  35.  
  36. ObjPtr isoClass;
  37.  
  38. /*Specifier for a part of the hex, which can be linked into a list*/
  39. typedef struct part 
  40.     {
  41.     struct part *next;    /*Next part in the list*/
  42.     int coords[3];        /*Coordinates of the part.  
  43.                   For a vertex, all three must be 0 or 1
  44.                   For an edge, exactly one must be 2
  45.                   For a face, two must be 2*/
  46.     } Part, *PartPtr;
  47.  
  48. /*Parts for using*/
  49. #define NPARTS    1580
  50. Part parts[NPARTS];
  51. PartPtr freePart, tempFree;
  52.  
  53. #define MAXNPOLYS 4
  54.  
  55. PartPtr presetPolys[256][MAXNPOLYS];
  56.  
  57. /*Macros for allocating parts*/
  58. #define NEWPART    (freePart ?                               \
  59.          (tempFree = freePart, freePart = freePart -> next, tempFree) \
  60.          : newp(Part))
  61.  
  62. #define FREEPART(p) p -> next = freePart; freePart = p;
  63.  
  64.  
  65. #ifdef PROTO
  66. Bool SectPart(PartPtr, PartPtr, PartPtr);
  67. Bool IsoHex(real fieldHex[2][2][2], real surfVal, ObjPtr pic,
  68.         Bool lowInside);
  69. #else
  70. Bool SectPart();
  71. Bool IsoHex();
  72. #endif
  73.  
  74. /*Macro to insert a node after one node in the linked list*/
  75. #define INSERTAFTER(list, node)                        \
  76.     {                                \
  77.         node -> next = list -> next;                 \
  78.         list -> next = node;                    \
  79.     }
  80.  
  81. /*Macro to dispose of a node after the current node*/
  82. #define DISPOSEAFTER(list)                        \
  83.     {                                \
  84.         PartPtr next;                        \
  85.         next = list -> next -> next;                \
  86.         FREEPART(list -> next);                    \
  87.         list -> next = next;                    \
  88.     }
  89.  
  90. /*Values for vertexState*/
  91. #define WITHIN    0        /*Within the isosurface*/
  92. #define WITHOUT 1        /*Outside the isosurface*/
  93. #define USED    2        /*Already used in the polygon*/
  94.  
  95. Bool EdgeOnFace(e, f)
  96. PartPtr e, f;
  97. /*Returns true if edge e is on face f*/
  98. {
  99.     register int k;
  100.     for (k = 0; k < 3; ++k)
  101.     {
  102.     if (e -> coords[k] == 2 && f -> coords[k] != 2)
  103.     {
  104.         return false;
  105.     }
  106.     }
  107.     return true;
  108. }
  109.  
  110. Bool SectPart(s1, s2, d)
  111. PartPtr s1, s2, d;
  112. /*Returns true if s1 intersects s2 and, if so, puts the intersection in d
  113.   The intersection of two edges is a vertex; the intersection of two faces
  114.   is an edge*/
  115. {
  116.     register int k;
  117.   
  118.     for (k = 0; k < 3; ++k)
  119.     {
  120.     if (s1 -> coords[k] == 2)
  121.     {
  122.         d -> coords[k] = s2 -> coords[k];
  123.     }
  124.     else if (s2 -> coords[k] == 2)
  125.     {
  126.         d -> coords[k] = s1 -> coords[k];
  127.     }
  128.     else if (s1 -> coords[k] != s2 -> coords[k])
  129.     {
  130.         return false;
  131.     }
  132.     else
  133.     {
  134.         d -> coords[k] = s1 -> coords[k];
  135.     }
  136.     }
  137.     return true;
  138. }
  139.  
  140. void UnionPart(s1, s2, d)
  141. PartPtr s1, s2, d;
  142. /*Puts in d the union of parts s1 and s2*/
  143. {
  144.     register int k;
  145.  
  146.     for (k = 0; k < 3; ++k)
  147.     {
  148.     if (s1 -> coords[k] == s2 -> coords[k])
  149.     {
  150.         d -> coords[k] = s1 -> coords[k];
  151.     }
  152.     else
  153.     {
  154.         d -> coords[k] = 2;
  155.     }
  156.     }
  157. }
  158.  
  159. Bool ThirdEdge(s1, s2, d)
  160. PartPtr s1, s2, d;
  161. /*Given 2 edges, s1, and s2 which meet at a point, forms the third edge d*/
  162. {
  163.     register int k;
  164.  
  165.     for (k = 0; k < 3; ++k)
  166.     {
  167.     if (s1 -> coords[k] == 2)
  168.     {
  169.         d -> coords[k] = s2 -> coords[k];
  170.     }
  171.     else if (s2 -> coords[k] == 2)
  172.     {
  173.         d -> coords[k] = s1 -> coords[k];
  174.     }
  175.     else if (s1 -> coords[k] != s2 -> coords[k])
  176.     {
  177.         return false;
  178.     }
  179.     else
  180.     {
  181.         d -> coords[k] = 2;
  182.     }
  183.     }
  184.     return true;
  185. }
  186.  
  187. void Other2Edges(e, v, e1, e2)
  188. PartPtr e, v, e1, e2;
  189. /*Given an edge e and a vertex v on that edge, puts in e1 and e2 the other
  190.   two edges which intersect that vertex*/
  191. {
  192.     register int index;                /*Index into which coordinate*/
  193.     
  194.     for (index = 0; index < 3; ++index)
  195.     {
  196.     if (e -> coords[index] != 2)
  197.     {
  198.         e1 -> coords[0] = v -> coords[0];
  199.         e1 -> coords[1] = v -> coords[1];
  200.         e1 -> coords[2] = v -> coords[2];
  201.         e1 -> coords[index] = 2;
  202.         break;
  203.     }
  204.     }
  205.  
  206.     for (index = 2; index >= 0; --index)
  207.     {
  208.     if (e -> coords[index] != 2)
  209.     {
  210.         e2 -> coords[0] = v -> coords[0];
  211.         e2 -> coords[1] = v -> coords[1];
  212.         e2 -> coords[2] = v -> coords[2];
  213.         e2 -> coords[index] = 2;
  214.         break;
  215.     }
  216.     }
  217. }
  218.  
  219. #ifdef PROTO
  220. Bool CacheIso(int initState[2][2][2], int whichPreset, Bool lowInside)
  221. #else
  222. Bool CacheIso(initState, whichPreset, lowInside)
  223. int initState[2][2][2];
  224. int whichPreset;
  225. Bool lowInside;
  226. #endif
  227. /*Calculates a series of polygons that make an isosurface through fieldHex
  228.   at surfVal, storing it in whichPreset.  If lowInside, assumes
  229.   that what is lower than surfVal is inside, otherwise what is higher*/
  230. {
  231.     int i, j, k;            /*Counters*/
  232.     int vertexState[2][2][2];        /*State of each vertex in the hex*/
  233.     int verticesLeft;            /*Number of within vertices left*/
  234.     Bool retVal = true;            /*Return value*/
  235.     int whichPolygon;            /*Which polygon of this polygon cache*/
  236.  
  237.     /*Figure out what the state of the vertices is*/
  238.     verticesLeft = 0;
  239.     for (i = 0; i < 2; ++i)
  240.     {
  241.     for (j = 0; j < 2; ++j)
  242.     {
  243.         for(k = 0; k < 2; ++k)
  244.         {
  245.         if (initState[i][j][k])
  246.         {
  247.             vertexState[i][j][k] = lowInside ? WITHIN : WITHOUT;
  248.         }
  249.         else
  250.         {
  251.             vertexState[i][j][k] = lowInside ? WITHOUT : WITHIN;
  252.         }
  253.         if (vertexState[i][j][k] == WITHIN)
  254.         {
  255.             ++verticesLeft;
  256.         }
  257.         }
  258.     }
  259.     }
  260.  
  261.     /*If all vertices are within, don't bother trying*/
  262.     if (verticesLeft == 8)
  263.     {
  264.     return true;
  265.     }
  266.  
  267.     whichPolygon = 0;
  268.     /*Create isosurfaces until no more are needed*/
  269.     while (verticesLeft)
  270.     {
  271.     int i0, j0, k0;            /*Index of starting vertex*/
  272.     PartPtr polygon;        /*The polygon around the isosurface,
  273.                       circular linked list*/
  274.     register PartPtr oneEdge;    /*Temporary edge*/
  275.  
  276.     /*Pick an initial vertex around which to build the surface*/
  277.     for (i = 0; i < 2; ++i)
  278.     {
  279.         for (j = 0; j < 2; ++j)
  280.         {
  281.         for (k = 0; k < 2; ++k)
  282.         {
  283.             if (vertexState[i][j][k] == WITHIN)
  284.             {
  285.             i0 = i;
  286.             j0 = j;
  287.             k0 = k;
  288.             }
  289.         }
  290.         }
  291.     }
  292.  
  293.     /*Now, i0, j0, and k0 are set to the initial vertex*/
  294.  
  295.     /*Create the first polygon and account for that vertex*/
  296.     oneEdge = NEWPART;
  297.     oneEdge -> next = oneEdge;
  298.     oneEdge -> coords[0] = 2;
  299.     oneEdge -> coords[1] = j0;
  300.     oneEdge -> coords[2] = k0;
  301.     polygon = oneEdge;
  302.  
  303.     /*Use checkerboard rule to determine initial order*/
  304.     if (i0 ^ j0 ^ k0)
  305.     {
  306.         oneEdge = NEWPART;
  307.         oneEdge -> coords[0] = i0;
  308.         oneEdge -> coords[1] = j0;
  309.         oneEdge -> coords[2] = 2;
  310.         INSERTAFTER(polygon, oneEdge);
  311.  
  312.         oneEdge = NEWPART;
  313.         oneEdge -> coords[0] = i0;
  314.         oneEdge -> coords[1] = 2;
  315.         oneEdge -> coords[2] = k0;
  316.         INSERTAFTER(polygon, oneEdge);
  317.     }
  318.     else
  319.     {
  320.         oneEdge = NEWPART;
  321.         oneEdge -> coords[0] = i0;
  322.         oneEdge -> coords[1] = 2;
  323.         oneEdge -> coords[2] = k0;
  324.         INSERTAFTER(polygon, oneEdge);
  325.  
  326.         oneEdge = NEWPART;
  327.         oneEdge -> coords[0] = i0;
  328.         oneEdge -> coords[1] = j0;
  329.         oneEdge -> coords[2] = 2;
  330.         INSERTAFTER(polygon, oneEdge);
  331.     }
  332.  
  333.     vertexState[i0][j0][k0] = USED;
  334.     --verticesLeft;
  335.  
  336.     /*Now go through and expand the polygon*/
  337.     for (;;)
  338.     {
  339.         PartPtr runner;        /*Runner through the polygon*/
  340.  
  341.         /*First see if there is any potential vertex, common to two 
  342.           existing adjacent edges which can combine them into one*/
  343.         runner = polygon;
  344.         do
  345.         {
  346.         Part intersection;    /*Intersection between two edges*/
  347.         if (SectPart(runner, runner -> next, &intersection))
  348.         {
  349.             /*intersection contains a vertex where two edges meet*/
  350.             if (vertexState[intersection . coords[0]]
  351.                    [intersection . coords[1]]
  352.                    [intersection . coords[2]] == WITHIN)
  353.             {
  354.             /*It's a valid candidate; gobble it up*/
  355.  
  356.             ThirdEdge(runner, runner -> next, runner);
  357.             if (runner -> next == polygon)
  358.             {
  359.                 polygon = runner;
  360.             }
  361.             DISPOSEAFTER(runner);
  362.             vertexState[intersection . coords[0]]
  363.                                    [intersection . coords[1]]
  364.                                    [intersection . coords[2]] = USED;
  365.             --verticesLeft;
  366.             goto tryAgain;
  367.             }
  368.         }
  369.         runner = runner -> next;
  370.         } while (runner != polygon);
  371.  
  372.         /*Now see if an edge can be stretched over a neighboring vertex*/
  373.         runner = polygon;
  374.         do
  375.         {
  376.         int index;        /*Index into which coord is Free*/
  377.         Part vertex;
  378.  
  379.         /*Start off with vertex = edge*/
  380.         vertex . coords[0] = runner -> coords[0];
  381.         vertex . coords[1] = runner -> coords[1];
  382.         vertex . coords[2] = runner -> coords[2];
  383.         
  384.         /*Figure out which coordinate is Free*/
  385.         for (index = 0; index < 3; ++index)
  386.         {
  387.             if (runner -> coords[index] == 2)
  388.             {
  389.             break;
  390.             }
  391.         }
  392.  
  393.         /*Now try both vertices that share that edge*/
  394.         for (vertex . coords[index] = 0;
  395.              vertex . coords[index] < 2;
  396.              ++(vertex . coords[index]))
  397.         {
  398.  
  399.             if (vertexState[vertex . coords[0]]
  400.                    [vertex . coords[1]]
  401.                    [vertex . coords[2]] == WITHIN)
  402.             {
  403.             /*Yes, it's good!  Snap it over it.*/
  404.             Part lastFace;        /*Second face to climb over*/
  405.             Part temp;        /*Temporary part*/
  406.             Part edge1, edge2;    /*Candidate edges*/
  407.  
  408.  
  409.             /*Determine the last face the new line will climb over
  410.               to determine the order of edges*/
  411.             UnionPart(runner, runner -> next, &lastFace);
  412.  
  413.             /*Start with the edge that does not intersect lastFace*/
  414.             Other2Edges(runner, &vertex, &edge1, &edge2);
  415.             
  416.             oneEdge = NEWPART;
  417.             if (EdgeOnFace(&edge1, &lastFace))
  418.             {
  419.                 /*Start with edge2*/
  420.                 runner -> coords[0] = edge2 . coords[0];
  421.                 runner -> coords[1] = edge2 . coords[1];
  422.                 runner -> coords[2] = edge2 . coords[2];
  423.  
  424.                 oneEdge -> coords[0] = edge1 . coords[0];
  425.                 oneEdge -> coords[1] = edge1 . coords[1];
  426.                 oneEdge -> coords[2] = edge1 . coords[2];
  427.             }
  428.             else
  429.             {
  430.                 /*Start with edge1*/
  431.                 runner -> coords[0] = edge1 . coords[0];
  432.                 runner -> coords[1] = edge1 . coords[1];
  433.                 runner -> coords[2] = edge1 . coords[2];
  434.  
  435.                 oneEdge -> coords[0] = edge2 . coords[0];
  436.                 oneEdge -> coords[1] = edge2 . coords[1];
  437.                 oneEdge -> coords[2] = edge2 . coords[2];
  438.             }
  439.  
  440.             INSERTAFTER(runner, oneEdge);
  441.  
  442.             vertexState[vertex . coords[0]]
  443.                                    [vertex . coords[1]]
  444.                                    [vertex . coords[2]] = USED;
  445.             --verticesLeft;
  446.             goto tryAgain;
  447.             }
  448.         }
  449.  
  450.         runner = runner -> next;
  451.         } while (runner != polygon);
  452.  
  453. #if 0
  454.         /*Failed the two tests.  See if there are two identical intersections between two used
  455.         vertices to slide*/
  456.         runner = polygon;
  457.         do
  458.         {
  459.         int index;        /*Index into which coord is Free*/
  460.         Part vertex;
  461.  
  462.         /*Start off with vertex = edge*/
  463.         vertex . coords[0] = runner -> coords[0];
  464.         vertex . coords[1] = runner -> coords[1];
  465.         vertex . coords[2] = runner -> coords[2];
  466.         
  467.         /*Figure out which coordinate is Free*/
  468.         for (index = 0; index < 3; ++index)
  469.         {
  470.             if (runner -> coords[index] == 2)
  471.             {
  472.             break;
  473.             }
  474.         }
  475.  
  476.         /*Try both endpoints to see if they're both USED*/
  477.         vertex . coords[index] = 0;
  478.         if (vertexState[vertex . coords[0]]
  479.                                [vertex . coords[1]]
  480.                                [vertex . coords[2]] == USED)
  481.         {
  482.             /*Looks good so far*/
  483.             vertex . coords[index] = 1;
  484.             if (vertexState[vertex . coords[0]]
  485.                                    [vertex . coords[1]]
  486.                                    [vertex . coords[2]] == USED)
  487.             {
  488.             /*Excellent! Found one.  Now search for the other one, if any*/
  489.             PartPtr other;
  490.             other = runner -> next;
  491.             do
  492.             {
  493.                 if (runner -> coords[0] == other -> coords[0] &&
  494.                 runner -> coords[1] == other -> coords[1] &&
  495.                 runner -> coords[2] == other -> coords[2])
  496.                 {
  497.                 /*Hot diggity sclotos!*/
  498.                 break;
  499.                 }
  500.                 other = other -> next;
  501.             } while (other != runner);
  502.             if (other != runner)
  503.             {
  504.                 /*We found a mate!  Slide the two*/
  505.  
  506.                 runner -> coords[0] = other -> next -> coords[0];
  507.                 runner -> coords[1] = other -> next -> coords[1];
  508.                 runner -> coords[2] = other -> next -> coords[2];
  509.  
  510.                 other -> coords[0] = runner -> next -> coords[0];
  511.                 other -> coords[1] = runner -> next -> coords[1];
  512.                 other -> coords[2] = runner -> next -> coords[2];
  513.                 break;
  514.             }
  515.             }
  516.         }
  517.  
  518.         runner = runner -> next;
  519.         } while (runner != polygon);
  520. #endif
  521.         break;
  522.         
  523. tryAgain:   ;
  524.     }
  525.  
  526.     if (polygon)
  527.     {
  528.         PartPtr curNode;        /*The current node in the polygon*/
  529.         PartPtr next;        /*The next node in the polygon*/
  530.         int whichVertex;        /*The current vertex*/
  531.         int k;
  532.  
  533.         whichVertex = 0;
  534.  
  535.         /*Test the polygon for sanity*/
  536.         curNode = polygon;
  537.         do
  538.         {
  539.             register int n;        /*Dimension counters*/
  540.         register int ic, jc, kc;
  541.  
  542.         int t1, t2;
  543.  
  544.         /*Determine along which side this intersection lies*/
  545.         ic = curNode -> coords[0];
  546.         jc = curNode -> coords[1];
  547.         kc = curNode -> coords[2];
  548.         if (ic == 2)
  549.         {
  550.             /*i is Free*/
  551.             t1 = initState[0][jc][kc];
  552.             t2 = initState[1][jc][kc];
  553.             if (t1 == t2)
  554.             {
  555.             retVal = false;
  556.             goto zapPolygon;
  557.             }
  558.         }
  559.         else if (jc == 2)
  560.         {
  561.             /*j is Free*/
  562.             t1 = initState[ic][0][kc];
  563.             t2 = initState[ic][1][kc];
  564.             if (t1 == t2)
  565.             {
  566.             retVal = false;
  567.             goto zapPolygon;
  568.             }
  569.         }
  570.         else if (kc == 2)
  571.         {
  572.             /*k is Free*/
  573.             t1 = initState[ic][jc][0];
  574.             t2 = initState[ic][jc][1];
  575.             if (t1 == t2)
  576.             {
  577.             retVal = false;
  578.             goto zapPolygon;
  579.             }
  580.         }
  581.         else
  582.         {
  583.             ReportError("CacheIso", "Error, neither ic nor kc nor jc Free");
  584.         }
  585.  
  586.         curNode = curNode -> next;
  587.         } while (curNode != polygon);
  588.  
  589.         if (whichPolygon >= MAXNPOLYS) ReportError("CacheIso", "Too many polys");
  590.         presetPolys[whichPreset][whichPolygon] = polygon;
  591.         ++whichPolygon;
  592.         continue;
  593. zapPolygon:
  594.         /*Get rid of the polygon*/
  595.         curNode = polygon;
  596.         next = curNode -> next;
  597.         do
  598.         {
  599.         curNode = next;
  600.         next = curNode -> next;
  601.         FREEPART(curNode);
  602.         } while (curNode != polygon);
  603.     }
  604.     }
  605.     return retVal;
  606. }
  607.  
  608. #ifdef PROTO
  609. Bool NewIsoHex(real *bottomField, real *topField,
  610.     VertexPtr *bottomVertices, VertexPtr *topVertices,
  611.     long index[], long dims[],
  612.     real surfVal, ObjPtr pic,
  613.     Bool lowInside)
  614. #else
  615. Bool NewIsoHex(bottomField, topField,
  616.     bottomVertices, topVertices,
  617.     index, dims,
  618.     surfVal, pic,
  619.     lowInside)
  620. real *bottomField; real *topField;
  621. VertexPtr *bottomVertices; VertexPtr *topVertices;
  622. long index[]; long dims[];
  623. real surfVal; ObjPtr pic;
  624. Bool lowInside;
  625. #endif
  626. /*Calculates a series of polygons that make an isosurface.*/
  627. {
  628.     int i, j, t;                /*Counters*/
  629.     real sample;
  630.     int whichCase;
  631.     register long offset;
  632.  
  633.     /*Calculate the case*/
  634.     whichCase = 0;
  635.     for (i = 0; i < 2; ++i)
  636.     {
  637.     for (j = 0; j < 2; ++j)
  638.     {
  639.         offset = (index[0] + i) * dims[1] + (index[1] + j);
  640.  
  641.         whichCase = whichCase << 1;
  642.  
  643.         sample = bottomField[offset];
  644.         if (sample == missingData) return true;
  645.         if (lowInside ? (sample <= surfVal) : (sample >= surfVal))
  646.         {
  647.         whichCase += 1;
  648.         }
  649.  
  650.         whichCase = whichCase << 1;
  651.  
  652.         sample = topField[offset];
  653.         if (sample == missingData) return true;
  654.         if (lowInside ? (sample <= surfVal) : (sample >= surfVal))
  655.         {
  656.         whichCase += 1;
  657.         }
  658.     }
  659.     }
  660.  
  661.     i = index[0];
  662.     j = index[1];
  663.     for (t = 0; presetPolys[whichCase][t]; ++t)
  664.     {
  665.     PartPtr polygon;
  666.     polygon = presetPolys[whichCase][t];
  667.     if (polygon)
  668.     {
  669.         PartPtr curNode;        /*The current node in the polygon*/
  670.         PartPtr next;        /*The next node in the polygon*/
  671.         VertexPtr vertices[20];    /*The vertices into the polygon*/
  672.         int whichVertex;        /*Which vertex of the poly looking at*/
  673.  
  674.         whichVertex = 0;
  675.  
  676.         /*Emit the polygon*/
  677.         curNode = polygon;
  678.         do
  679.         {
  680.         register int ic, jc, kc;
  681.  
  682.         /*Determine along which side this intersection lies*/
  683.         ic = curNode -> coords[0];
  684.         jc = curNode -> coords[1];
  685.         kc = curNode -> coords[2];
  686.         if (ic == 2)
  687.         {
  688.             /*i is Free*/
  689.             offset = (i * dims[1] + (j + jc)) * 3;
  690.         }
  691.         else if (jc == 2)
  692.         {
  693.             /*j is Free*/
  694.             offset = ((i + ic) * dims[1] + j) * 3 + 1;
  695.         }
  696.         else if (kc == 2)
  697.         {
  698.             /*k is Free*/
  699.             offset = ((i + ic) * dims[1] + (j + jc)) * 3 + 2;
  700.         }
  701.         else
  702.         {
  703.             ReportError("NewIsoHex", "Error, neither ic nor kc nor jc Free");
  704.             return false;
  705.         }
  706.  
  707.         vertices[whichVertex] = (kc == 1) ? topVertices[offset] : bottomVertices[offset];
  708.  
  709.         if (vertices[whichVertex] == 0)
  710.         {
  711.             ReportError("NewIsoHex", "Vertex error");
  712.             return false;
  713.         }
  714.         ++whichVertex;
  715.  
  716.         curNode = curNode -> next;
  717.         } while (curNode != polygon);
  718.  
  719.         /*Now the polygon has been assembled; spit it out.*/
  720.         TesselateSPolyToPicture(pic, whichVertex, vertices);
  721.     }
  722.     }
  723. }
  724.  
  725. #define CALCDELTA(field, hc, pic, mic, pjc, mjc, pkc, mkc)
  726.  
  727. #undef    NEXTFIELD
  728. #undef    PREVFIELD
  729. #undef    PLUSIC
  730. #undef    MINUSIC
  731. #undef    PLUSJC
  732. #undef    MINUSJC
  733. #undef    PLUSKC
  734. #undef    MINUSKC
  735.  
  736.  
  737. #ifdef PROTO
  738. static void StuffIJKVertices(VertexPtr *elements, real *field, real *nextField, real deltas[], int fieldNum, int coord, long dims[3], long kVal, ObjPtr pic, real isoVal, Bool encloseLow)
  739. #else
  740. static void StuffIJKVertices(elements, field, nextField, deltas, fieldNum, coord, dims, kVal, pic, isoVal, encloseLow)
  741. VertexPtr *elements;
  742. real *field;
  743. int fieldNum;
  744. real *nextField;
  745. real deltas[];
  746. int coord;
  747. long dims[3];
  748. long kVal;
  749. ObjPtr pic;
  750. real isoVal;
  751. Bool encloseLow;
  752. #endif
  753. /*Makes a group of IJK vertices for field at k = kVal, sticking them in pic
  754.   Coord is the coordinates.  Stuffs the results into array.  If deltas, uses
  755.   deltas as a holding place for the deltas to make normals*/
  756. {
  757.     long index[3];
  758.     register long nElements;
  759.     register long i, j;
  760.     register long offset;
  761.     register real curVal, nextVal;
  762.     register real r;
  763.  
  764.     /*Empty the array*/
  765.     nElements = dims[0] * dims[1] * 3;
  766.     for (offset = 0; offset < nElements; ++offset)
  767.     {
  768.     elements[offset] = 0;
  769.     }
  770.  
  771.     /*Make I vertices*/
  772.     for (index[1] = 0; index[1] < dims[1]; ++index[1])
  773.     {
  774.     index[0] = 0;
  775.     index[2] = kVal;
  776.  
  777.     offset = index[1];
  778.  
  779.     curVal = field[offset];
  780.  
  781.     for (index[0] = 1; index[0] < dims[0]; ++index[0])
  782.     {
  783.         nextVal = field[offset + dims[1]];
  784.  
  785.         if (curVal != missingData && nextVal != missingData)
  786.         {
  787.         /*Maybe make some sort of a vertex*/
  788.  
  789.         if (BETWEENP(isoVal, curVal, nextVal))
  790.         {
  791.             /*Yes!  Vertex here!*/
  792.             register VertexPtr v;
  793.             register real location, location1;
  794.             real curCoord[3], nextCoord[3];
  795.  
  796.             v = NewVertex(pic, 0);
  797.             location = (isoVal - curVal) / (nextVal - curVal);
  798.             location1 = 1.0 - location;
  799.  
  800.             /*Get coordinates*/
  801.             --index[0];
  802.             curCoord[0] = SelectFieldComponent(coord, 0, index);
  803.             curCoord[1] = SelectFieldComponent(coord, 1, index);
  804.             curCoord[2] = SelectFieldComponent(coord, 2, index);
  805.             ++index[0];
  806.             nextCoord[0] = SelectFieldComponent(coord, 0, index);
  807.             nextCoord[1] = SelectFieldComponent(coord, 1, index);
  808.             nextCoord[2] = SelectFieldComponent(coord, 2, index);
  809.  
  810.             /*Set position*/
  811.             v -> position[0] = nextCoord[0] * location +
  812.                        curCoord[0] * location1;
  813.             v -> position[1] = nextCoord[1] * location +
  814.                        curCoord[1] * location1;
  815.             v -> position[2] = nextCoord[2] * location +
  816.                        curCoord[2] * location1;
  817.  
  818.             if (deltas)
  819.             {
  820.             /*Calculate normal from delta*/
  821.             --index[0];
  822. #define FIELDHERE    field
  823. #define NEXTFIELD    nextField
  824. #define HEREC        curCoord
  825. #define HEREV        curVal
  826. #define PLUSIC        nextCoord
  827. #include "ScianIsoCalcDelta.h"
  828. #undef FIELDHERE
  829. #undef NEXTFIELD
  830. #undef HEREC
  831. #undef HEREV
  832. #undef PLUSIC
  833.             ++index[0];
  834.  
  835. #define FIELDHERE    field
  836. #define NEXTFIELD    nextField
  837. #define HEREC        nextCoord
  838. #define HEREV        nextVal
  839. #define MINUSIC        curCoord
  840. #include "ScianIsoCalcDelta.h"
  841. #undef FIELDHERE
  842. #undef NEXTFIELD
  843. #undef HEREC
  844. #undef HEREV
  845. #undef MINUSIC
  846.  
  847.             /*Now that we have two deltas, make the normal*/
  848.             offset += dims[1];
  849.  
  850.             v -> normal[0] = location * deltas[offset * 3];
  851.             v -> normal[1] = location * deltas[offset * 3 + 1];
  852.             v -> normal[2] = location * deltas[offset * 3 + 2];
  853.  
  854.             offset -= dims[1];
  855.  
  856.             v -> normal[0] += location1 * deltas[offset * 3];
  857.             v -> normal[1] += location1 * deltas[offset * 3 + 1];
  858.             v -> normal[2] += location1 * deltas[offset * 3 + 2];
  859.  
  860.             r = sqrt(SQUARE(v -> normal[0]) +
  861.                  SQUARE(v -> normal[1]) +
  862.                  SQUARE(v -> normal[2]));
  863.             if (r > 0.0)
  864.             {
  865.                 r = 1.0 / r;
  866.                 v -> normal[0] *= r;
  867.                 v -> normal[1] *= r;
  868.                 v -> normal[2] *= r;
  869.  
  870.                 if (!encloseLow)
  871.                 {
  872.                 v -> normal[0] = -(v -> normal[0]);
  873.                 v -> normal[1] = -(v -> normal[1]);
  874.                 v -> normal[2] = -(v -> normal[2]);
  875.                 }
  876.             }
  877.             else
  878.             {
  879.                 /*Just give it plus z normal*/
  880.  
  881.                 v -> normal[0] = 0.0;
  882.                 v -> normal[1] = 0.0;
  883.                 v -> normal[2] = 1.0;
  884.             }
  885.             }
  886.             else
  887.             {
  888.             /*Just give it plus z normal*/
  889.  
  890.             v -> normal[0] = 0.0;
  891.             v -> normal[1] = 0.0;
  892.             v -> normal[2] = 1.0;
  893.             }
  894.  
  895.             elements[offset * 3] = v;
  896.         }
  897.         }
  898.         curVal = nextVal;
  899.         offset += dims[1];
  900.     }
  901.     }
  902.  
  903.     /*Make J vertices*/
  904.     for (index[0] = 0; index[0] < dims[0]; ++index[0])
  905.     {
  906.     index[1] = 0;
  907.     index[2] = kVal;
  908.  
  909.     offset = index[0] * dims[1];
  910.  
  911.     curVal = field[offset];
  912.  
  913.     for (index[1] = 1; index[1] < dims[1]; ++index[1])
  914.     {
  915.         nextVal = field[offset + 1];
  916.  
  917.         if (curVal != missingData && nextVal != missingData)
  918.         {
  919.         /*Maybe make some sort of a vertex*/
  920.  
  921.         if (BETWEENP(isoVal, curVal, nextVal))
  922.         {
  923.             /*Yes!  Vertex here!*/
  924.             register VertexPtr v;
  925.             register real location, location1;
  926.             real curCoord[3], nextCoord[3];
  927.  
  928.             v = NewVertex(pic, 0);
  929.             location = (isoVal - curVal) / (nextVal - curVal);
  930.             location1 = 1.0 - location;
  931.  
  932.             /*Get coordinates*/
  933.             --index[1];
  934.             curCoord[0] = SelectFieldComponent(coord, 0, index);
  935.             curCoord[1] = SelectFieldComponent(coord, 1, index);
  936.             curCoord[2] = SelectFieldComponent(coord, 2, index);
  937.             ++index[1];
  938.             nextCoord[0] = SelectFieldComponent(coord, 0, index);
  939.             nextCoord[1] = SelectFieldComponent(coord, 1, index);
  940.             nextCoord[2] = SelectFieldComponent(coord, 2, index);
  941.  
  942.             /*Set position*/
  943.             v -> position[0] = nextCoord[0] * location +
  944.                        curCoord[0] * location1;
  945.             v -> position[1] = nextCoord[1] * location +
  946.                        curCoord[1] * location1;
  947.             v -> position[2] = nextCoord[2] * location +
  948.                        curCoord[2] * location1;
  949.  
  950.             if (deltas)
  951.             {
  952.             /*Calculate normal from delta*/
  953.             --index[1];
  954. #define FIELDHERE    field
  955. #define NEXTFIELD    nextField
  956. #define HEREC        curCoord
  957. #define HEREV        curVal
  958. #define PLUSJC        nextCoord
  959. #include "ScianIsoCalcDelta.h"
  960. #undef FIELDHERE
  961. #undef NEXTFIELD
  962. #undef HEREC
  963. #undef HEREV
  964. #undef PLUSJC
  965.             ++index[1];
  966.  
  967. #define FIELDHERE    field
  968. #define NEXTFIELD    nextField
  969. #define HEREC        nextCoord
  970. #define HEREV        nextVal
  971. #define MINUSJC        curCoord
  972. #include "ScianIsoCalcDelta.h"
  973. #undef FIELDHERE
  974. #undef NEXTFIELD
  975. #undef HEREC
  976. #undef HEREV
  977. #undef MINUSJC
  978.  
  979.             /*Now that we have two deltas, make the normal*/
  980.             ++offset;
  981.  
  982.             v -> normal[0] = location * deltas[offset * 3];
  983.             v -> normal[1] = location * deltas[offset * 3 + 1];
  984.             v -> normal[2] = location * deltas[offset * 3 + 2];
  985.  
  986.             --offset;
  987.  
  988.             v -> normal[0] += location1 * deltas[offset * 3];
  989.             v -> normal[1] += location1 * deltas[offset * 3 + 1];
  990.             v -> normal[2] += location1 * deltas[offset * 3 + 2];
  991.  
  992.             r = sqrt(SQUARE(v -> normal[0]) +
  993.                  SQUARE(v -> normal[1]) +
  994.                  SQUARE(v -> normal[2]));
  995.             if (r > 0.0)
  996.             {
  997.                 r = 1.0 / r;
  998.                 v -> normal[0] *= r;
  999.                 v -> normal[1] *= r;
  1000.                 v -> normal[2] *= r;
  1001.  
  1002.                 if (!encloseLow)
  1003.                 {
  1004.                 v -> normal[0] = -(v -> normal[0]);
  1005.                 v -> normal[1] = -(v -> normal[1]);
  1006.                 v -> normal[2] = -(v -> normal[2]);
  1007.                 }
  1008.             }
  1009.             else
  1010.             {
  1011.                 /*Just give it plus z normal*/
  1012.  
  1013.                 v -> normal[0] = 0.0;
  1014.                 v -> normal[1] = 0.0;
  1015.                 v -> normal[2] = 1.0;
  1016.             }
  1017.             }
  1018.             else
  1019.             {
  1020.             /*Just give it plus z normal*/
  1021.  
  1022.             v -> normal[0] = 0.0;
  1023.             v -> normal[1] = 0.0;
  1024.             v -> normal[2] = 1.0;
  1025.             }
  1026.  
  1027.             elements[offset * 3 + 1] = v;
  1028.         }
  1029.         }
  1030.         curVal = nextVal;
  1031.         ++offset;
  1032.     }
  1033.     }
  1034.  
  1035.     /*Make K vertices*/
  1036.     index[0] = 0;
  1037.     index[1] = 0;
  1038.     index[2] = kVal;
  1039.  
  1040.     if (kVal < dims[2])
  1041.     for (index[0] = 0; index[0] < dims[0]; ++index[0])
  1042.     {
  1043.     offset = index[0] * dims[1];
  1044.     for (index[1] =0; index[1] < dims[1]; ++index[1])
  1045.     {
  1046.         index[2] = kVal;
  1047.         curVal = field[offset];
  1048.  
  1049.         index[2] = kVal + 1;
  1050.         nextVal = nextField[offset];
  1051.  
  1052.         if ((curVal != missingData) &&
  1053.         (nextVal != missingData) &&
  1054.         BETWEENP(isoVal, curVal, nextVal))
  1055.         {
  1056.             /*Yes!  Vertex here!*/
  1057.             register VertexPtr v;
  1058.             register real location, location1;
  1059.             real curCoord[3], nextCoord[3];
  1060.  
  1061.             v = NewVertex(pic, 0);
  1062.             location = (isoVal - curVal) / (nextVal - curVal);
  1063.             location1 = 1.0 - location;
  1064.  
  1065.             /*Get coordinates*/
  1066.             --index[2];
  1067.             curCoord[0] = SelectFieldComponent(coord, 0, index);
  1068.             curCoord[1] = SelectFieldComponent(coord, 1, index);
  1069.             curCoord[2] = SelectFieldComponent(coord, 2, index);
  1070.  
  1071.             ++index[2];
  1072.             nextCoord[0] = SelectFieldComponent(coord, 0, index);
  1073.             nextCoord[1] = SelectFieldComponent(coord, 1, index);
  1074.             nextCoord[2] = SelectFieldComponent(coord, 2, index);
  1075.  
  1076.             /*Set position*/
  1077.             v -> position[0] = nextCoord[0] * location +
  1078.                        curCoord[0] * location1;
  1079.             v -> position[1] = nextCoord[1] * location +
  1080.                        curCoord[1] * location1;
  1081.             v -> position[2] = nextCoord[2] * location +
  1082.                        curCoord[2] * location1;
  1083.  
  1084.             if (deltas)
  1085.             {
  1086.             /*Calculate normal from delta*/
  1087.             --index[2];
  1088. #define FIELDHERE    field
  1089. #define NEXTFIELD    nextField
  1090. #define HEREC        curCoord
  1091. #define PLUSKC        nextCoord
  1092. #define HEREV        curVal
  1093. #include "ScianIsoCalcDelta.h"
  1094. #undef FIELDHERE
  1095. #undef NEXTFIELD
  1096. #undef HEREC
  1097. #undef PLUSKC
  1098. #undef HEREV
  1099.             v -> normal[0] = location1 * deltas[offset * 3];
  1100.             v -> normal[1] = location1 * deltas[offset * 3 + 1];
  1101.             v -> normal[2] = location1 * deltas[offset * 3 + 2];
  1102.  
  1103.             deltas[offset * 3] = missingData;
  1104.             deltas[offset * 3 + 1] = missingData;
  1105.             deltas[offset * 3 + 2] = missingData;
  1106.  
  1107.             ++index[2];
  1108.  
  1109. #define FIELDHERE    nextField
  1110. #define PREVFIELD    field
  1111. #define HEREC        nextCoord
  1112. #define MINUSKC        curCoord
  1113. #define HEREV        nextVal
  1114. #include "ScianIsoCalcDelta.h"
  1115. #undef FIELDHERE
  1116. #undef PREVFIELD
  1117. #undef HEREC
  1118. #undef MINUSKC
  1119. #undef HEREV
  1120.             /*Now that we have two deltas, make the normal*/
  1121.  
  1122.             v -> normal[0] += location * deltas[offset * 3];
  1123.             v -> normal[1] += location * deltas[offset * 3 + 1];
  1124.             v -> normal[2] += location * deltas[offset * 3 + 2];
  1125.  
  1126.             r = sqrt(SQUARE(v -> normal[0]) +
  1127.                  SQUARE(v -> normal[1]) +
  1128.                  SQUARE(v -> normal[2]));
  1129.             if (r > 0.0)
  1130.             {
  1131.                 r = 1.0 / r;
  1132.                 v -> normal[0] *= r;
  1133.                 v -> normal[1] *= r;
  1134.                 v -> normal[2] *= r;
  1135.  
  1136.                 if (!encloseLow)
  1137.                 {
  1138.                 v -> normal[0] = -(v -> normal[0]);
  1139.                 v -> normal[1] = -(v -> normal[1]);
  1140.                 v -> normal[2] = -(v -> normal[2]);
  1141.                 }
  1142.             }
  1143.             else
  1144.             {
  1145.                 /*Just give it plus z normal*/
  1146.  
  1147.                 v -> normal[0] = 0.0;
  1148.                 v -> normal[1] = 0.0;
  1149.                 v -> normal[2] = 1.0;
  1150.             }
  1151.             }
  1152.             else
  1153.             {
  1154.             /*Just give it plus z normal*/
  1155.  
  1156.             v -> normal[0] = 0.0;
  1157.             v -> normal[1] = 0.0;
  1158.             v -> normal[2] = 1.0;
  1159.             }
  1160.  
  1161.             elements[offset * 3 + 2] = v;
  1162.         }
  1163.         else
  1164.         {
  1165.             /*Haven't made deltas for next time, must kill them*/
  1166.             if (deltas)
  1167.             {
  1168.                 deltas[offset * 3] = missingData;
  1169.                 deltas[offset * 3 + 1] = missingData;
  1170.                 deltas[offset * 3 + 2] = missingData;
  1171.             }
  1172.         }
  1173.  
  1174.         curVal = nextVal;
  1175.         ++offset;
  1176.     }
  1177.     }
  1178. }
  1179.  
  1180. static ObjPtr SetIsoMainDataset(visObj, dataSet)
  1181. ObjPtr visObj, dataSet;
  1182. /*Sets field to be the main field in visObj*/
  1183. {
  1184.    SetVar(visObj, MAINDATASET, dataSet);
  1185.    return ObjTrue;
  1186. }
  1187.  
  1188. static ObjPtr AddIsoControls(object, panelContents)
  1189. ObjPtr object, panelContents;
  1190. /*Adds controls appropriate to an isosurface to panelContents*/
  1191. {
  1192.     int left, right, top;
  1193.     real initValue;        /*The initial value of a slider*/
  1194.     ObjPtr var;            /*A variable of some sort*/
  1195.     ObjPtr slider;        /*A slider*/
  1196.     ObjPtr textBox;        /*A text box*/
  1197.     ObjPtr corral;        /*An icon corral*/
  1198.     int width;            /*Width of the area we can use*/
  1199.     ObjPtr name;        /*The name of the field*/
  1200.     ObjPtr icon;        /*An icon that represents the field*/
  1201.     ObjPtr isoField;        /*The field of the isosurface.*/
  1202.     ObjPtr defaultIcon;        /*The default icon of the field*/
  1203.     real db[2];            /*Bounds of the isosurface*/
  1204.     ObjPtr minMax;        /*Array containing min and max*/
  1205.     ObjPtr checkBox;        /*Temporary check box*/
  1206.     ObjPtr titleBox;        /*Temporary title box*/
  1207.     ObjPtr radioGroup;        /*Group of radio buttons*/
  1208.     ObjPtr mainDataset;        /*The main dataset of the field*/
  1209.     ObjPtr scale;        /*Scale for slider*/
  1210.  
  1211.     width = CWINWIDTH - 2 * CORRALBORDER - CWINCORRALWIDTH;
  1212.     left = MAJORBORDER;
  1213.     top = CWINHEIGHT - MAJORBORDER;
  1214.  
  1215.     /*Put in the isosurface corral at the left*/
  1216.     corral = NewIconCorral(NULLOBJ,
  1217.                left, left + ONECORRALWIDTH,
  1218.                top - ONECORRALHEIGHT,
  1219.                top,
  1220.                0);
  1221.     SetVar(corral, SINGLECORRAL, ObjTrue);
  1222.     SetVar(corral, TOPDOWN, ObjTrue);
  1223.     SetVar(corral, NAME, NewString("Source Field"));
  1224.     PrefixList(panelContents, corral);
  1225.     SetVar(corral, HELPSTRING,
  1226.     NewString("This corral shows the dataset that is being used to make \
  1227. the isosurface.  To replace it with another dataset, drag the icon of the other \
  1228. dataset into this corral."));
  1229.     SetVar(corral, PARENT, panelContents);
  1230.     SetVar(corral, REPOBJ, object);
  1231.     SetMethod(corral, DROPINCONTENTS, DropInMainDatasetCorral);
  1232.  
  1233.     /*Create the iso source text box*/
  1234.     textBox = NewTextBox(left, left + ONECORRALWIDTH, 
  1235.              top - ONECORRALHEIGHT - TEXTBOXSEP - TEXTBOXHEIGHT,
  1236.              top - ONECORRALHEIGHT - TEXTBOXSEP,
  1237.              0, "Isosurface Field Text", "Source Field");
  1238.     PrefixList(panelContents, textBox);
  1239.     SetVar(textBox, PARENT, panelContents);
  1240.     SetTextAlign(textBox, CENTERALIGN);
  1241.     top = top - ONECORRALHEIGHT - TEXTBOXSEP - TEXTBOXHEIGHT - MAJORBORDER;
  1242.     right = left + ONECORRALWIDTH;
  1243.  
  1244.     /*Put in an icon that represents the field*/
  1245.     isoField = GetObjectVar("AddIsoControls", object, MAINDATASET);
  1246.     if (!isoField) return ObjFalse;
  1247.     MakeVar(isoField, MINMAX);
  1248.     minMax = GetVar(isoField, MINMAX);
  1249.     if (!minMax) return ObjFalse;
  1250.  
  1251.     Array2CArray(db, minMax);
  1252.     while (mainDataset = GetVar(isoField, MAINDATASET))
  1253.     {
  1254.     isoField = mainDataset;
  1255.     }
  1256.  
  1257.     name = GetVar(isoField, NAME);
  1258.     defaultIcon = GetVar(isoField, DEFAULTICON);
  1259.     if (defaultIcon)
  1260.     {
  1261.     icon = NewObject(defaultIcon, 0);
  1262.     SetVar(icon, NAME, name);
  1263.     SetVar(icon, REPOBJ, isoField);
  1264.     }
  1265.     else
  1266.     {
  1267.     icon = NewIcon(0, 0, ICONQUESTION, GetString(name));
  1268.     }
  1269.     SetVar(icon, ICONLOC, NULLOBJ);
  1270.     DropIconInCorral(corral, icon);
  1271.  
  1272.     ChooseGoodRange(&(db[0]), &(db[1]));
  1273.  
  1274.     /*Create the isovalue slider*/
  1275.     slider = TemplateSlider(IsosurfaceTemplate, "Isovalue", SCALE);
  1276.     if (!slider)
  1277.     {
  1278.     return ObjFalse;
  1279.     }
  1280.     PrefixList(panelContents, slider);
  1281.     SetVar(slider, PARENT, panelContents);
  1282.     SetVar(slider, HELPSTRING,
  1283.     NewString("This slider controls the isosurface value shown by an \
  1284. isosurface visualization object.  Move the indicator to a new desired isosurface value and \
  1285. and a new isosurface will be calculated.  If you don't want the isosurface \
  1286. to be recalculated, you can turn off the icon representing the isosurface \
  1287. before you do this.  You can get multiple isosurfaces by duplicating the \
  1288. isosurface icon."));
  1289.  
  1290.     scale = TemplateScale(IsosurfaceTemplate, "Iso scale",
  1291.             SO_LEFT, false);
  1292.     PrefixList(panelContents, scale);
  1293.     SetVar(scale, PARENT, panelContents);
  1294.     LinkScale(scale, slider);
  1295.     
  1296.     MakeVar(object, ISOVAL);
  1297.     var = GetRealVar("AddIsoControls", object, ISOVAL);
  1298.     if (var)
  1299.     {
  1300.     initValue = GetReal(var);
  1301.     }
  1302.     else
  1303.     {
  1304.     initValue = (db[0] + db[1]) * 0.5;
  1305.     SetVar(object, ISOVAL, NewReal(initValue));
  1306.     }
  1307.     SetSliderRange(slider, db[1], db[0], 0.0);
  1308.  
  1309.     AssocDirectControlWithVar(slider, object, ISOVAL);
  1310.     SetTrackNot(slider, true);
  1311.  
  1312.     /*Create the value legend box*/
  1313.     textBox = TemplateTextBox(IsosurfaceTemplate, "Isovalue Text", ONE_LINE, "Isovalue:");
  1314.     PrefixList(panelContents, textBox);
  1315.     SetVar(textBox, PARENT, panelContents);
  1316.  
  1317.     top -= TEXTBOXHEIGHT + MINORBORDER;
  1318.  
  1319.     /*Create the value text box*/
  1320.     textBox = TemplateTextBox(IsosurfaceTemplate, "Isovalue Readout", EDITABLE + WITH_PIT + ONE_LINE, "Isovalue");
  1321.     PrefixList(panelContents, textBox);
  1322.     SetVar(textBox, PARENT, panelContents);
  1323.     SetTextAlign(textBox, RIGHTALIGN);
  1324.     SliderReadout(slider, textBox);
  1325.  
  1326.     /*Make radio buttons for enclosing*/
  1327.  
  1328.     radioGroup = NewRadioButtonGroup("Enclose Data");
  1329.     SetVar(radioGroup, HELPSTRING,
  1330.     NewString("These radio buttons allow you to specify whether high or \
  1331. are to be enclosed within the isosurfaces.  This information is used to determine \
  1332. the direction of surface normals and is also used to improve the behavior of the \
  1333. isosurface routine in ambiguous cases.")); 
  1334.  
  1335.     checkBox = TemplateRadioButton(IsosurfaceTemplate, "High Values");
  1336.     SetVar(checkBox, HELPSTRING,
  1337.     NewString("This button specifies that the shapes of the isosurface  \
  1338. enclose values higher than the isosurface value."));
  1339.     AddRadioButton(radioGroup, checkBox);
  1340.  
  1341.     checkBox = TemplateRadioButton(IsosurfaceTemplate, "Low Values");
  1342.     SetVar(checkBox, HELPSTRING,
  1343.     NewString("This button specifies that the shapes of the isosurface \
  1344. enclose values lower than the isosurface value."));
  1345.     AddRadioButton(radioGroup, checkBox);
  1346.  
  1347.     /*Title box around it*/
  1348.     titleBox = TemplateTitleBox(IsosurfaceTemplate, "Enclose");
  1349.     PrefixList(panelContents, titleBox);
  1350.     SetVar(titleBox, PARENT, panelContents);
  1351.  
  1352.     /*Add the radio button group*/
  1353.     PrefixList(panelContents, radioGroup);
  1354.     SetVar(radioGroup, PARENT, panelContents);
  1355.  
  1356.     /*Set its value*/
  1357.     if (!GetIntVar("AddIsoControls", object, ENCLOSELOW))
  1358.     {
  1359.     SetVar(object, ENCLOSELOW, NewInt(0));
  1360.     }
  1361.     AssocDirectControlWithVar(radioGroup, object, ENCLOSELOW);
  1362.  
  1363.     /*Get Normals*/
  1364.     titleBox = TemplateTitleBox(IsosurfaceTemplate, "Get Normals");
  1365.     PrefixList(panelContents, titleBox);
  1366.     SetVar(titleBox, PARENT, panelContents);
  1367.  
  1368.     radioGroup = NewRadioButtonGroup("Normal Radio");
  1369.     SetVar(radioGroup, HELPSTRING,
  1370.     NewString("These radio buttons allow you to specify whether the \
  1371. vertex normals of the resulting surface are taken from the geometry or \
  1372. from the data.  Taking the normals from the data is more accurate and \
  1373. produces a smoother surface.  \
  1374. Taking the normals from the geometry is faster.")); 
  1375.  
  1376.     checkBox = TemplateRadioButton(IsosurfaceTemplate, "From Data");
  1377.     SetVar(checkBox, HELPSTRING,
  1378.     NewString("This button specifies that the surface normals of the \
  1379. isosurface be taken from the gradient of the data, run through trilinear \
  1380. interpolation."));
  1381.     AddRadioButton(radioGroup, checkBox);
  1382.  
  1383.     checkBox = TemplateRadioButton(IsosurfaceTemplate, "From Geometry");
  1384.     SetVar(checkBox, HELPSTRING,
  1385.     NewString("This button specifies that the surface normals of the \
  1386. isosurface be calculated from the geometry using Phong lighting."));
  1387.     AddRadioButton(radioGroup, checkBox);
  1388.  
  1389.     /*Add the radio button group*/
  1390.     PrefixList(panelContents, radioGroup);
  1391.     SetVar(radioGroup, PARENT, panelContents);
  1392.  
  1393.     AssocDirectControlWithVar(radioGroup, object, NORMALSFROM);
  1394.  
  1395.     return ObjTrue;
  1396. }
  1397.  
  1398. #define NHIST    20
  1399. #define HISTFRACT 0.25
  1400. #define MAXDESIRABLE    2000
  1401.  
  1402. static ObjPtr MakeIsoVal(object)
  1403. ObjPtr object;
  1404. /*Makes the isoval of an object*/
  1405. {
  1406.     real db[2];
  1407.     real isoVal;
  1408.     ObjPtr repObj;
  1409.     ObjPtr minMax;
  1410.     long histogram[NHIST];
  1411.     int k;
  1412.     int nTraversalDims;
  1413.     long *traversalDims;
  1414.     long *traversalSteps;
  1415.     long *index;
  1416.     real d;
  1417.     long total;
  1418.     long runningTotal;
  1419.     int subDivision;
  1420.     int advantage;
  1421.     int whichDim;
  1422.  
  1423.     MakeVar(object, MAINDATASET);
  1424.     repObj = GetObjectVar("MakeIsoVal", object, MAINDATASET);
  1425.     if (!repObj)
  1426.     {
  1427.     return ObjFalse;
  1428.     }
  1429.  
  1430.     MakeVar(repObj, MINMAX);
  1431.  
  1432.     minMax = GetFixedArrayVar("MakeIsoVal", repObj, MINMAX, 1, 2L);
  1433.  
  1434.     if (!minMax)
  1435.     {
  1436.     return ObjFalse;
  1437.     }
  1438.     Array2CArray(db, minMax);
  1439.     ChooseGoodRange(&(db[0]), &(db[1]));
  1440.  
  1441.     /*Now we have the range, get the iso value*/
  1442.     SetCurField(FIELD1, repObj);
  1443.  
  1444.     /*Get the information on traversing the dataset*/
  1445.     nTraversalDims = CountTraversalDims(FIELD1);
  1446.     if (nTraversalDims)
  1447.     {
  1448.     traversalDims = (long *) Alloc(sizeof(long) * nTraversalDims);
  1449.     traversalSteps = (long *) Alloc(sizeof(long) * nTraversalDims);
  1450.     index = (long *) Alloc(sizeof(long) * nTraversalDims);
  1451.     GetTraversalDims(FIELD1, traversalDims);
  1452.  
  1453.     /*Figure out how coarsely to traverse dataset*/
  1454.     runningTotal = 0;
  1455.     advantage = 1;
  1456.     for (k = 0; k < nTraversalDims; ++k)
  1457.     {
  1458.         runningTotal *= traversalDims[k];
  1459.         advantage *= 2;
  1460.     }
  1461.     subDivision = 1;
  1462.     while (runningTotal > MAXDESIRABLE)
  1463.     {
  1464.         ++subDivision;
  1465.         runningTotal /= advantage;
  1466.     }
  1467.  
  1468.     /*Zero the histogram*/
  1469.     for (k = 0; k < NHIST; ++k)
  1470.     {
  1471.         histogram[k] = 0;
  1472.     }
  1473.  
  1474.     /*Zero the index*/
  1475.     for (k = 0; k < nTraversalDims; ++k)
  1476.     {
  1477.         index[k] = 0;
  1478.     }
  1479.  
  1480.     /*Fill up the histogram*/
  1481.     do
  1482.     {
  1483.         d = SelectFieldScalar(FIELD1, index);
  1484.         if (d != missingData)
  1485.         {
  1486.         k = ((real) NHIST) * (d - db[0]) / (db[1] - db[0]);
  1487.         if (k < 0) k = 0;
  1488.         if (k >= NHIST) k = NHIST - 1;
  1489.         ++histogram[k];
  1490.         }
  1491.  
  1492.         /*Advance to next sample*/
  1493.         for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
  1494.         {
  1495.         if (traversalDims[whichDim] > 0)
  1496.         {
  1497.             index[whichDim] += subDivision;
  1498.             if (index[whichDim] >= traversalDims[whichDim])
  1499.             {
  1500.             index[whichDim] = 0;
  1501.             }
  1502.             else
  1503.             {
  1504.             break;
  1505.             }
  1506.         }
  1507.         }
  1508.     } while (whichDim < nTraversalDims); /*Break is based on advance*/
  1509.  
  1510.     /*Now choose from histogram*/
  1511.     total = 0;
  1512.     for (k = 0; k < NHIST; ++k)
  1513.     {
  1514.         total += histogram[k];
  1515.     }
  1516.     runningTotal = 0;
  1517.     for (k = NHIST - 1; k > 1; --k)
  1518.     {
  1519.         runningTotal += histogram[k];
  1520.         if (runningTotal > HISTFRACT * total)
  1521.         {
  1522.         break;
  1523.         }
  1524.     }
  1525.     ++k;
  1526.     if (k > NHIST - 1) --k;
  1527.  
  1528.     isoVal = k * (db[1] - db[0]) / NHIST + db[0];
  1529.  
  1530.     SAFEFREE(index);
  1531.     SAFEFREE(traversalDims);
  1532.     SAFEFREE(traversalSteps);
  1533.     }
  1534.     else
  1535.     {
  1536.     isoVal = (db[0] + db[1]) * 0.5;
  1537.     }
  1538.  
  1539.     SetVar(object, ISOVAL, NewReal(isoVal));
  1540.     return ObjTrue;
  1541. }
  1542.  
  1543. Bool Escape()
  1544. {
  1545.     int x, y;
  1546.     if (runningScript)
  1547.     {
  1548.     return false;
  1549.     }
  1550.     else
  1551.     {
  1552.     return Mouse(&x, &y);
  1553.     }
  1554. }
  1555.  
  1556. ObjPtr NewMakeIsoSurface(object)
  1557. ObjPtr object;
  1558. /*Makes a SURFACE in isosurface vis object*/
  1559. {
  1560.     ObjPtr repObj;
  1561.     long dims[3];
  1562.     long index[3];
  1563.     double time;
  1564.     ObjPtr var;
  1565.     real isoVal;
  1566.     Bool encloseLow;
  1567.     ObjPtr surface;
  1568.     int isoLevel;
  1569.     ObjPtr lastIJK, nextIJK;
  1570.     ObjPtr lastData, nextData, nextPlusData;
  1571.     Bool useDeltas;
  1572.     ObjPtr deltas;
  1573.  
  1574.     time = Clock();
  1575.  
  1576.     repObj = GetObjectVar("MakeIsoSurface", object, MAINDATASET);
  1577.     if (!repObj)
  1578.     {
  1579.     return ObjFalse;
  1580.     }
  1581.  
  1582.     /*Get the iso value*/
  1583.     MakeVar(object, ISOVAL);
  1584.     var = GetRealVar("MakeIsoSurface", object, ISOVAL);
  1585.     if (!var)
  1586.     {
  1587.     return ObjFalse;
  1588.     }
  1589.     isoVal = GetReal(var);
  1590.  
  1591.     encloseLow = GetPredicate(object, ENCLOSELOW);
  1592.  
  1593.     useDeltas = GetPredicate(object, NORMALSFROM) ? false : true;
  1594.  
  1595.     LongOperation();
  1596.  
  1597.     /*Get the isolevel to start at*/
  1598.     var = GetVar(object, ISOLEVEL);
  1599.     if (var)
  1600.     {
  1601.     isoLevel = GetInt(var);
  1602.     }
  1603.     else
  1604.     {
  1605.     isoLevel = 0;
  1606.     }
  1607.  
  1608.     SetCurField(FIELD1, repObj);
  1609.     SetCurForm(FORMFIELD, repObj);
  1610.     if (CountTraversalDims(FIELD1) != 3)
  1611.     {
  1612.     ReportError("MakeIsoSurface", "Wrong number of traversal dimensions");
  1613.     return ObjFalse;
  1614.     }
  1615.     GetTraversalDims(FIELD1, dims);
  1616.  
  1617.     if (isoLevel)
  1618.     {
  1619.     /*This is just a continuation of a previous incarnation*/
  1620.     if (isoLevel >= dims[2])
  1621.     {
  1622.         /*Escape clause; just return*/
  1623.         return ObjTrue;
  1624.     }
  1625.     surface = GetVar(object, SURFACE);
  1626.     if (!surface)
  1627.     {
  1628.         surface = NewPicture();
  1629.         SetVar(object, SURFACE, surface);
  1630.     }
  1631.  
  1632.     lastIJK = GetVar(object, LASTIJKVERTICES);
  1633.     if (!lastIJK)
  1634.     {
  1635.         ReportError("MakeIsoSurface", "No LASTIJKVERTICES");
  1636.         return ObjFalse;
  1637.     }
  1638.  
  1639.     nextIJK = GetVar(object, NEXTIJKVERTICES);
  1640.     if (!nextIJK)
  1641.     {
  1642.         ReportError("MakeIsoSurface", "No NEXTIJKVERTICES");
  1643.         return ObjFalse;
  1644.     }
  1645.  
  1646.     lastData = GetVar(object, LASTDATA);
  1647.     if (!lastData)
  1648.     {
  1649.         ReportError("MakeIsoSurface", "No LASTDATA");
  1650.         return ObjFalse;
  1651.     }
  1652.  
  1653.     nextData = GetVar(object, NEXTDATA);
  1654.     if (!nextData)
  1655.     {
  1656.         ReportError("MakeIsoSurface", "No NEXTDATA");
  1657.         return ObjFalse;
  1658.     }
  1659.  
  1660.     deltas = GetVar(object, DELTAS);
  1661.     }
  1662.     else
  1663.     {
  1664.     long arrayDims[3];
  1665.     real *elements;
  1666.     long nElements;
  1667.  
  1668.     surface = NewPicture();
  1669.     SetVar(surface, REPOBJ, object);
  1670.     SetVar(object, SURFACE, surface);
  1671.  
  1672.     arrayDims[0] = dims[0];
  1673.     arrayDims[1] = dims[1];
  1674.     arrayDims[2] = 3;
  1675.  
  1676.     /*Make data arrays*/
  1677.     lastData = NewArray(AT_REAL, 2, arrayDims);
  1678.     SetVar(object, LASTDATA, lastData);
  1679.  
  1680.     nextData = NewArray(AT_REAL, 2, arrayDims);
  1681.     SetVar(object, NEXTDATA, nextData);
  1682.  
  1683.     if (useDeltas)
  1684.     {
  1685.         deltas = NewArray(AT_REAL, 3, arrayDims);
  1686.  
  1687.         /*Make them all missing*/
  1688.         nElements = arrayDims[0] * arrayDims[1] * arrayDims[2];
  1689.         elements = ELEMENTS(deltas);
  1690.         do
  1691.         {
  1692.         --nElements;
  1693.         elements[nElements] = missingData;
  1694.         } while (nElements);
  1695.     }
  1696.     else
  1697.     {
  1698.         deltas = NULLOBJ;
  1699.     }
  1700.  
  1701.     /*Fill in the data*/
  1702.     index[0] = 0;
  1703.     index[1] = 0;
  1704.     index[2] = 0;
  1705.     elements = ELEMENTS(lastData);
  1706.     StuffIJPlane(elements, FIELD1, 0, index);
  1707.  
  1708.     /*Get last data*/
  1709.     index[0] = 0;
  1710.     index[1] = 0;
  1711.     index[2] = 0;
  1712.     StuffIJPlane((real *) ELEMENTS(lastData), FIELD1, 0, index);
  1713.  
  1714.     /*Get next data*/
  1715.     index[0] = 0;
  1716.     index[1] = 0;
  1717.     index[2] = 1;
  1718.     StuffIJPlane((real *) ELEMENTS(nextData), FIELD1, 0, index);
  1719.  
  1720.     /*Make IJK arrays*/
  1721.     lastIJK = NewArray(AT_POINTER, 3, arrayDims);
  1722.     SetVar(object, LASTIJKVERTICES, lastIJK);
  1723.  
  1724.     nextIJK = NewArray(AT_POINTER, 3, arrayDims);
  1725.     SetVar(object, NEXTIJKVERTICES, nextIJK);
  1726.  
  1727.     /*Fill in just the lastIJK*/
  1728.     StuffIJKVertices(ELEMENTS(lastIJK), ELEMENTS(lastData), ELEMENTS(nextData), deltas ? ELEMENTS(deltas) : 0, FIELD1, FORMFIELD, dims, 0, surface, isoVal, encloseLow);
  1729.     }
  1730.  
  1731.     nextPlusData = NewArray(AT_REAL, 2, dims);
  1732.  
  1733.     /*Go through and do the isosurface*/
  1734.     while (isoLevel < dims[2] - 1)
  1735.     {
  1736.     index[0] = 0;
  1737.     index[1] = 0;
  1738.     index[2] = isoLevel + 2;
  1739.  
  1740.     /*Get next plus data*/
  1741.     if (isoLevel < dims[2] - 2)
  1742.     {
  1743.         StuffIJPlane((real *) ELEMENTS(nextPlusData), FIELD1, 0, index);
  1744.     }
  1745.  
  1746.     /*Get next IJK's*/
  1747.     StuffIJKVertices(ELEMENTS(nextIJK), ELEMENTS(nextData), ELEMENTS(nextPlusData), deltas ? ELEMENTS(deltas) : 0, FIELD1, FORMFIELD, dims, isoLevel + 1, surface, isoVal, encloseLow);
  1748.  
  1749.     /*Now go through and make the isosurface*/
  1750.     for (index[0] = 0; index[0] < dims[0] - 1; ++index[0])
  1751.     {
  1752.         for (index[1] = 0; index[1] < dims[1] - 1; ++index[1])
  1753.         {
  1754.         NewIsoHex(ELEMENTS(lastData), ELEMENTS(nextData),
  1755.             ELEMENTS(lastIJK), ELEMENTS(nextIJK),
  1756.             index, dims,
  1757.             isoVal, surface, encloseLow);
  1758.         }
  1759.     }
  1760.  
  1761.     /*Shift down ijk vertices*/
  1762.     var = lastIJK;
  1763.     lastIJK = nextIJK;
  1764.     nextIJK = var;
  1765.  
  1766.     /*Shift down data*/
  1767.     var = lastData;
  1768.     lastData = nextData;
  1769.     nextData = nextPlusData;
  1770.     nextPlusData = var;
  1771.  
  1772.     /*Go to next isoLevel*/
  1773.     ++isoLevel;
  1774.  
  1775. #if 0
  1776.     /*Insert escape code here*/
  1777.     if (Escape())
  1778.     {
  1779.         SetVar(object, LASTIJKVERTICES, lastIJK);
  1780.         SetVar(object, NEXTIJKVERTICES, nextIJK);
  1781.         SetVar(object, LASTDATA, lastData);
  1782.         SetVar(object, NEXTDATA, nextData);
  1783.         SetVar(object, ISOLEVEL, NewInt(isoLevel));
  1784.         SetVar(object, DELTAS, deltas);
  1785.         DeferMessage(object, MAKEISOSURFACE);
  1786.         return ObjFalse;
  1787.     }
  1788. #endif
  1789.     }
  1790.  
  1791.     if (!deltas)
  1792.     {
  1793.     CalcPictureNormals(surface);
  1794.     }
  1795.  
  1796.     SetVar(object, SURFACE, surface);
  1797.     MakeVar(object, CPALETTE);
  1798.     SetVar(surface, CPALETTE, GetVar(object, CPALETTE));
  1799.  
  1800.     return ObjTrue;
  1801. }
  1802.  
  1803. static ObjPtr MakeIsoSurfaceLater(object)
  1804. ObjPtr object;
  1805. /*Sets up to make an isosurface later*/
  1806. {
  1807.     FuncTyp method;
  1808.  
  1809.     SetVar(object, ISOLEVEL, NewInt(0));
  1810.     method = GetMethod(object, MAKEISOSURFACE);
  1811.  
  1812.     if (method)
  1813.     {
  1814.     (*method)(object);
  1815.     }
  1816.     return ObjTrue;
  1817. }
  1818.  
  1819. static ObjPtr IsoInit(iso)
  1820. ObjPtr iso;
  1821. /*Initializes an isosurface*/
  1822. {
  1823.     ObjPtr mainDataset;
  1824.     long nDims;
  1825.     long *dims;
  1826.     ObjPtr retVal;
  1827.  
  1828.     MakeVar(iso, MAINDATASET);
  1829.     mainDataset = GetVar(iso, MAINDATASET);
  1830.     if (!mainDataset) return NULLOBJ;
  1831.  
  1832.     SetCurForm(FORMFIELD, mainDataset);
  1833.     nDims = CountTraversalDims(FORMFIELD);
  1834.     if (nDims != 3) return NULLOBJ;
  1835.  
  1836.     dims = Alloc(3 * sizeof(long));
  1837.     GetTraversalDims(FORMFIELD, dims);
  1838.  
  1839.     retVal = ObjTrue;
  1840.  
  1841.     if (dims[0] <= 1) retVal = ObjFalse;
  1842.     if (dims[1] <= 1) retVal = ObjFalse;
  1843.     if (dims[2] <= 1) retVal = ObjFalse;
  1844.  
  1845.     return retVal;
  1846. }
  1847.  
  1848. void InitIsosurfaces()
  1849. {
  1850.     ObjPtr icon;
  1851.     int i, j, k, n;
  1852.     int initState[2][2][2];
  1853.     int *statePtr;
  1854.     int stateFlags, tempFlags;
  1855.  
  1856.     /*Initialize the Free space*/
  1857.     parts[0] . next = 0;
  1858.     for (k = 1; k < NPARTS; ++k)
  1859.     {
  1860.     parts[k] . next = &(parts[k - 1]);
  1861.     }
  1862.     freePart = &(parts[NPARTS - 1]);
  1863.  
  1864.     /*Zero the polygons*/
  1865.     for (i = 0; i < 256; ++i)
  1866.     {
  1867.     for (j = 0; j < MAXNPOLYS; ++j)
  1868.     {
  1869.         presetPolys[i][j] = 0;
  1870.     }
  1871.     }
  1872.  
  1873.     /*Fill the polygons*/
  1874.     for (stateFlags = 0; stateFlags < 256; ++stateFlags)
  1875.     {
  1876.     /*Fill initState based on stateFlag*/
  1877.     tempFlags = stateFlags;
  1878.     statePtr = &(initState[0][0][0]);
  1879.     for (k = 0; k < 8; ++k)
  1880.     {
  1881.         *statePtr++ = (tempFlags & 128) ? 1 : 0;
  1882.         tempFlags = tempFlags << 1;
  1883.     }
  1884.  
  1885.     if (!CacheIso(initState, stateFlags, false))
  1886.     {
  1887.         CacheIso(initState, stateFlags, true);
  1888.     }
  1889.     }
  1890.  
  1891.     /*Class for an isosurface*/
  1892.     isoClass = NewObject(visSurface, 0);
  1893.     AddToReferenceList(isoClass);
  1894.     SetVar(isoClass, NAME, NewString("Isosurface"));
  1895.     SetMethod(isoClass, ADDCONTROLS, AddIsoControls);
  1896.     SetVar(isoClass, DEFAULTICON, icon = NewObject(visIcon, 0));
  1897.     SetVar(isoClass, NORMALSFROM, NewInt(0));
  1898.     SetVar(icon, WHICHICON, NewInt(ICONISOSURFACE));
  1899.     SetVar(icon, NAME, NewString("Isosurface"));
  1900.     SetVar(icon, HELPSTRING,
  1901.     NewString("This icon represents an isosurface.  \
  1902. The isosurface object generates surfaces of constant value through 3-D \
  1903. scalar datasets defined over structured grids.  Some day it will work with \
  1904. nonstructured grids, as well."));
  1905.     icon = NewIcon(0, 0, ICONISOSURFACE, "Isosurface");
  1906.     SetVar(icon, HELPSTRING,
  1907.     NewString("Select this icon to show controls for the isosurface, such \
  1908. as the isosurface level and the treatment of missing data."));
  1909.     SetVar(isoClass, CONTROLICON, icon);
  1910.  
  1911.     SetMethod(isoClass, SETMAINDATASET, SetIsoMainDataset);
  1912.     SetMethod(isoClass, INITIALIZE, IsoInit);
  1913.     SetVar(isoClass, SAVEEXTENSION, NewString("isoSurface"));
  1914.     SetVar(isoClass, HANDLEMISSING, NewInt(0));
  1915.     SetMethod(isoClass, ISOVAL, MakeIsoVal);
  1916.     DeclareDependency(isoClass, SURFACE, ISOVAL);
  1917.     DeclareDependency(isoClass, SURFACE, ENCLOSELOW);
  1918.     SetVar(isoClass, ENCLOSELOW, NewInt(0));
  1919.     DeclareIndirectDependency(isoClass, SURFACE, MAINDATASET, CHANGED);
  1920.     DeclareDependency(isoClass, SURFACE, MAINDATASET);
  1921.     DeclareDependency(isoClass, SURFACE, NORMALSFROM);
  1922.     SetMethod(isoClass, MAKEISOSURFACE, NewMakeIsoSurface);
  1923.     SetMethod(isoClass, SURFACE, MakeIsoSurfaceLater);
  1924.  
  1925.     DefineVisMapping(DS_HASFORM | DS_HASFIELD, 3, 3, -1, isoClass);
  1926.     DefineVisMapping(DS_HASFORM | DS_HASFIELD | DS_VECTOR, 3, 3, -1, isoClass);
  1927. }
  1928.  
  1929. void KillIsosurfaces()
  1930. {
  1931.     DeleteThing(isoClass);
  1932. }
  1933.